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Abstract 

We present a bifurcation analysis of a normal form for travelling waves in one-dimensional 
excitable media. The normal form which has been recently proposed on phenomenological 
grounds is given in form of a differential delay equation. The normal form exhibits a symme- 
try preserving Hopf bifurcation which may coalesce with a saddle-node in a Bogdanov-Takens 
point, and a symmetry breaking spatially inhomogeneous pitchfork bifurcation. We study 
here the Hopf bifurcation for the propagation of a single pulse in a ring by means of a cen- 
ter manifold reduction, and for a wave train by means of a multiscale analysis leading to a 
real Ginzburg-Landau equation as the corresponding amplitude equation. Both, the center 
manifold reduction and the multiscale analysis show that the Hopf bifurcation is always sub- 
critical independent of the parameters. This may have links to cardiac alternans which have 
so far been believed to be stable oscillations emanating from a supercritical bifurcation. We 
discuss the implications for cardiac alternans and revisit the instability in some excitable 
media where the oscillations had been believed to be stable. In particular, we show that 
our condition for the onset of the Hopf bifurcation coincides with the well known restitution 
condition for cardiac alternans. 
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Excitable media are abundant in nature. Examples range from small scale systems 
such as intracellular calcium waves to large scale systems such as cardiac tissue. 
There exists a plethora of models describing excitable media, each of those par- 
ticular to the microscopic details of the underlying biological, chemical or physical 
system. However, excitable media have certain features which are common to all 
these systems. In a recent paper we introduced a normal form for travelling waves 
in one-dimensional excitable media which contains all bifurcations occurring in ex- 
citable media. The normal form consists of a delay-differential equation and is ap- 
plicable to systems which are close to the saddle-node bifurcation of the traveling 
wave. Although the normal form has so far only been proposed on phenomeno- 
logical grounds and has not yet been rigorously derived, its parameters could be 
successfully fitted to some real excitable media with good quantitative agreement. 
In this work we perform a bifurcation analysis of the occurring Hopf bifurcation. 
This may have important consequences for cardiac dynamics and the understanding 
of arrhythmias, in particular of alternans. Alternans describe the scenario in car- 
diac tissue whereby action potential durations are alternating periodically between 
short and long periods. There is an increased interest in alternans because they are 
believed to trigger spiral wave breakup in cardiac tissue and to be a precursor to 
ventricular fibrillation. So far these alternans were believed to be stable. However 
within the normal form we show that the Hopf bifurcation is actually subcritical 
suggesting that the resulting oscillations may be unstable. 

1 Introduction 

Examples of excitable media are frequently found in biological and chemical systems. Prominent 
examples are cardiac and neural tissue [U [2] , slime mold colonies in a starving environment [3] 
and intracellular calcium waves [1]. There are two defining features of excitable media which 
are crucial to enable effective signal transmission in biological systems such as cardiac or neural 
tissue: threshold behaviour and relaxation to a stable rest state. The threshold behaviour as- 
sures that only for large enough stimuli a signal is produced whereas small perturbation decay 
immediately. For super-threshold perturbations a signal will decay only after a long excursion 
- called action potential in the context of cardiac dynamics - back to its stable rest state. This 
relaxation allows for repeated stimulation which is essential for wave propagation in cardiac and 
neural tissue. 

Typical solutions in one-dimensional excitable media are wave trains. The wavelength L can 
range from L = oo corresponding to an isolated pulse to a minimal value Lc below which prop- 
agation fails. Besides wave trains rotating spiral waves may form in two dimensional excitable 
media, and scroll waves in three dimensional excitable media O [6l [71 [H [9] . In the context of 
cardiac excitable media propagation failure of these solutions is often linked to clinical situa- 
tions. In particular the break-up of spiral waves has been associated with pathological cardiac 
arrhythmias [lO]. Spiral waves may be created in cardiac tissue when wave trains propagate 
through inhomogeneities of the cardiac tissue. A reentrant spiral may move around an anatomi- 
cal obstacle or around a region of partially or totally inexcitable tissue. Once created they drive 
the heart at a rate much faster than normal sinus rhythm and cause tachycardia. If these spiral 
waves then subsequently breakup into multiple drifting and meandering spirals and disintegrate 
into a disorganized state, fibrillation may occur with a possible fatal result for the patient, es- 



1 



pecially when occurring in the ventricles. It is therefore of great interest to understand the 
transition from one reentrant spiral to the disorganized collection of complex reentrant path- 
ways. Rather than investigating the full two-dimensional problem of spiral wave break-up which 
would include interactions of numerous wave arms, one can study some aspects of spiral wave 
breakup by looking at a one-dimensional slice of a spiral i.e. at a one-dimensional wave train 
A pulse circulating around a one-dimensional ring constitutes the simplest model for a 
spiral rotating around an anatomical obstacle. Such models concentrate solely on the dynamics 
close to the anatomical obstacle and ignore the influence of the dynamics of the spiral arms. 

Experimentally this problem has been studied since the beginning of the last century. In [13] 
the circulation of an electrical pulse around a ring-shaped piece of atrial cardiac tissue from 
a tortoise heart was investigated as a model for reentrant activity. Recently, the experiments 
in [H] investigated the dynamics of pulse circulation around a ring-shaped piece of myocardial 
tissue from a dog heart. Remarkably, it was found that oscillations of the circulation period of 
the pulse occurred leading to conduction block and subsequent termination of reentry. 

This indicates that a one-dimensional oscillatory instability, named alternans, may be the mech- 
anism triggering spiral wave breakup \15\ \TT\ [T6\ \T7\ I18j . This instability occurs when the 
circulation time of the pulse around the ring is below a certain specific threshold. Below this 
threshold alternans arise and action potential durations are alternating periodically between 
short and long periods. We make here the distinction between alternans which are mediated 
by external periodic stimulations [Il[20l[2ll[22l[23l[23l[25l[26l[27! and alternans which have 
a purely dynamical origin. We call the latter ones dynamical alternans. It is these dynamical 
alternans with which this paper is concerned with. We will analyze dynamical alternans using 
a normal form proposed in [28j for travelling waves in one-dimensional excitable media. The 
stability of dynamical alternans will be determined by a center manifold theory and by multiple 
scale analysis. Typically oscillatory instabilities arise in one-dimensional media via Hopf bifur- 
cations. These Hopf bifurcations may be supercritical resulting in sustained stable oscillations or 
subcritical leading to a collapse of the oscillations and possibly of the pulse solution. Alternans 
are widely believed to originate via a supercritical bifurcation [HI [171 Ell [30]. This believe is 
based on numerical simulations of certain models for cardiac dynamics. However we note that 
the above mentioned experiments [T3| 114] show that the occurrence of oscillations leads to a 
subsequent termination of the pulse, which is not suggestive of a supercritical bifurcation. 

Our main result for a single pulse and for a wave train in a ring is that in the framework of 
our normal form alternans arise as a subcritical Hopf bifurcation. This result is in contrast to 
the common belief that alternans are stable oscillations. We corroborate our result by numeri- 
cal simulations of a model for cardiac tissue in which previous numerical simulations suggested 
the occurrence of stable oscillations. We will see that previous numerical experiments have not 
been performed sufficiently long to reveal the subcritical character of the Hopf bifurcation. The 
subcritical character is in agreement with the above mentioned experiments [131 [T4j and may 
explain why alternans often trigger spiral wave breakup and are associated with cardiac failures. 
The normal form describes excitable media in parameter regions where the system is close to the 
saddle-node of the travelling wave. This is the case for models such as the FitzHugh-Nagumo 
model [31] , and is typically the case when the Hopf bifurcation occurs close to the saddle-node 
bifurcation of the isolated pulse. In particular the normal form is valid for those classes of 
excitable media (or parameter regions of a particular excitable medium) where the activator 
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Figure 1: Plot of the activator u (continuous line) and the inhibitor v (dashed line) for the mod- 
ified Barkley system (|70p showing how the activator u weakly interacts with the exponentially 
decaying tail of the inhibitor v. The parameters are a = 0.22, Us = 0.1, e = 0.03755. The ring 
length is L = 245 which is close to the critical value Lc for the saddle-node bifurcation. 



weakly interacts with the preceding inhibitor which exponentially decays towards the homoge- 
neous rest state. However for other models such as the model by Echebarria and Karma [21] 
there exist parameter regions where upon decreasing the length of the ring the pulse solution 
is driven away from the solitary pulse solution and does not weakly interact with the inhibitor. 
Our normal form cannot describe these oscillations and we present numerical simulations where 
oscillations are indeed stable in Section 15.31 

Recently we have constructed a normal form for travelling waves in one-dimensional excitable 
media which takes the form of a delay differential equation [28] (see ([I]) and The construc- 
tion is based on the well-known observation that the interaction of a pulse with the inhibitor of 
the preceding pulse modifies the generic saddle-node bifurcation of an isolated pulse. In Fig. [1] 
we illustrate this scenario for a modified Barkley model [32|. The normal form exhibits a rich 
bifurcation behaviour which we could verify by numerically simulating partial differential equa- 
tion models of excitable media. Besides the well known saddle-node bifurcations for isolated 
pulses and for periodic wave trains the normal form also exhibits a Hopf bifurcation and a 
symmetry breaking, spatially inhomogeneous pitchfork bifurcation. Moreover, the normal form 
shows that the saddle-node and the Hopf bifurcation are an unfolding of a Bogdanov-Takens 
point as previously suggested in [331 El]- The Hopf bifurcation is found to occur before the 
saddle-node bifurcation for a single pulse in a ring. For a wave train consisting of several pulses 
in a ring, the Hopf- and the saddle-node bifurcations occur after the symmetry breaking pitch- 
fork bifurcation in which every second pulse dies. We could verify these scenarios in numerical 
simulations of a modified Barkley-model [32] and the FitzHugh-Nagumo equations [31]. The 
normal form provides a unified framework to study all possible bifurcations of travelling waves 
in one-dimensional excitable media. 

We were able to determine the parameters of the normal form from numerical simulations of 
the modified Barkley model [321 134] . Using these numerically determined parameters we showed 
excellent agreement between the normal form and the full partial differential equation. We 
quantitatively described the Hopf bifurcation and the inhomogeneous pitchfork bifurcation with 
the normal form. Moreover, we were able to quantify the Bogdanov-Takens bifurcation. 
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Whereas the subcritical character of the pitchfork bifurcation had been estabhshed in [28], a 
detailed analysis of the Hopf bifurcation was missing. For example, an interesting question is 
whether the Hopf bifurcation is sub-critical (as numerically observed for the parameters chosen 
in [28] for the modified Barkley model), or whether it is possible to observe sustained stable 
oscillations. This has important implications for cardiac alternans as described above. In this 
paper we will analyze the Hopf bifurcation of the normal form in detail. We derive a normal form 
for the Hopf bifurcation which allows us to determine the stability of the bifurcating solutions 
close to criticality. 

Before we embark on the analytical investigation of the Hopf bifurcation and the implications 
for cardiac dynamics, we state that all conclusions drawn obviously depend on the validity of 
the normal form. So far the proposed normal form which we briefly review in Section [2] has 
not been rigorously derived for any excitable medium. However, we point out that in [28] we 
have shown good quantitative agreement with some real excitable media. We will discuss the 
limitations of our approach in more detail in Section [3 

In Section [2] we recall the normal form and some of its properties. In Section [3] we perform a 
center manifold reduction of the normal form to describe the character of the Hopf bifurcation 
for a single pulse in a ring. In Section |4] we look at the case where a pseudo continuum of 
modes undergo a Hopf bifurcation and derive in a multiple scale analysis a Ginzburg-Landau 
equation which allows us to study the stability of the Hopf bifurcation of a wave train. The 
paper concludes with Section [5] where we present results from numerical simulations, discuss 
the implications of our analysis to cardiac dynamics and make connections to previous studies 
on alternans. In particular, we show that our condition for the onset of the Hopf bifurcation 
coincides with the well known restitution condition for cardiac alternans. 



2 A normal form 

In |28j we introduced a normal form for a single pulse on a periodic domain with length L 

dtX = - gX^ - /3(7 + X{t - r) + ^iX{t)) , (1) 

where 

/? = /3oexp(-Kr) , (2) 

for positive /?, k, 7 and 71. Here X{t) may be for example the difference of the amplitude or 
the velocity of a pulse to its respective value at the saddle-node. The terms proportional to (3 
incorporate finite domain effects associated with the activator of an excitable medium running 
into its own inhibitor with speed cq after the temporal delay r = (L — i^)/co where v is the finite 
width of the pulse. Note that for /? = (i.e. for the isolated pulse with r 00) we recover the 
generic saddle-node bifurcation which is well known for excitable media. Numerical simulations 
of excitable media show that the bifurcations of a single propagating pulse in a ring are different 
from the bifurcations of a wave train consisting of several distinct pulses. In the case of a wave 
train with finite wave length where a pulse may run into the inhibitor created by its preceding 
pulse, we showed that it was sufficient to consider two alternating populations of pulses X and 
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Y. We derived the following extension 

dtX = -^-gX^ -p{-f + Y{t-T)+-fiX{t)) 

dtV = -fi-gY^ -(]{j + X{t-T)+^iY{t)) . (3) 

To avoid confusion we state that we use the term normal form here in two different contexts. 
Equations ([1]) and ([3]) are coined "normal form" as they are an attempt to describe the behaviour 
of travelling waves in generic one-dimensional excitable media. However, these normal forms 
exhibit a rather rich bifurcation behaviour. We will focus in Section [3] on the solution behaviour 
close to criticality. In that context we will speak of "normal forms" in the sense of bifurcation 
theory. 



2.1 Properties of the normal form 
2.1.1 Saddle-node bifurcation 

Equation ([T|) supports the following stationary solutions 

Xi,2 = + 7i) ± VPH-i- + 11? - 45(/^ + m ■ (4) 

It is readily seen that the upper solution branch is stable whereas the lower one is unstable [28]. 
The two solutions coalesce in a saddle-node bifurcation with 

XsN = -^{l+li) at ^sN = ^{l + li? -^1 . (5) 

One expects the parameter /3, which describes the coupling to the inhibitor, to be small (see 
[28]). This implies ^SN < 0, indicating that the saddle- node of a single pulse or of a periodic 
wave train on a finite ring occurs at smaller values of the bifurcation parameter ^ than for the 
isolated pulse. Hence, the bifurcation is shifted to the left when compared to the isolated pulse 
(see Figure [2]) . 

Besides this stationary saddle node bifurcation the normal form ([1]) also contains a Hopf bifur- 
cation. 



2.1.2 Hopf bifurcation 

Linearization of the normal form around the homogeneous solution X with respect to small 
perturbations of the form 6X exp At with X = a + iuj yields 

X + 2gX + /371 + /3e-^" = . (6) 

Besides the stationary saddle-node bifurcation ([5]) with A = 0, the linearization ([6]) also reveals 
the existence of a Hopf bifurcation with X = itu. We readily find from ([6]) 

u) = /3sina;r (7) 
Xh = -|-(cosa;T + 7i) . (8) 
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Figure 2: Sketch of the bifurcation diagram for a single pulse in a ring showing a stationary 
saddle- node bifurcation (SN) and a subcritical Hopf bifurcation (HH). 



The first equation ^ allows us to formulate a necessary condition for the existence for a Hopf 
bifurcation 

/3r > 1 , 

i.e. if the coupling is strong enough and the pulse feels the presence of the inhibitor of the 
preceding pulse sufficiently strongly. The Hopf bifurcation occurs in parameter space before the 
saddle-node bifurcation and bifurcates from the upper stable branch Xi of (jl]) as is readily seen 
by observing Xh > XsN, independent of the value of /3. Hence we may equate Xh = Xi and 
solve for the bifurcation parameter. Setting fin = fisN — Sfi we find 

6fi = ^{l- cos ujt)^ . (9) 

The saddle-node bifurcation coalesces with the Hopf bifurcation in a codimension-2 Bogdanov- 
Takens point. At the Bogdanov-Takens point with fin = //5Ar we have wr = and uj = 0, i.e. 
the period of the oscillation goes to infinity. The amplitude at the Bogdanov-Takens point is 
readily determined by comparison of ([5]) with ([8]) for lot = 0. From ^ we infer that this occurs 
at pT = 1. We note that if /3r is large enough there can be arbitrary many solutions of d?]). 
We will discuss this scenario in Section HI 

In Figure [2] we show a schematic bifurcation diagram with the saddle-node bifurcation and 
the subcritical Hopf bifurcation for a single pulse in a ring. 

2.1.3 Spatially inhomogeneous pitchfork bifurcation 

When a group of several pulses in a ring is numerically simulated one observes that this wave train 
group does not undergo a symmetry preserving Hopf bifurcation on increasing the refractoriness, 
but instead develops a symmetry breaking, spatially inhomogeneous instability whereby every 
second pulse dies. 

The normal form is able to predicted and quantitatively describe this scenario [28] . The 
system ([3]) for wave trains supports two types of solutions. Besides the homogeneous solution 
dH), Xh = Yh = Xi, which may undergo a saddle- node bifurcation as described by ([5]), there 
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exists another stationary solution, an alternating mode Xa and 1^, with 

Xa = -Ya + ^{l-ll) . (10) 

9 

Associated with this alternating solution is a pitchfork bifurcation at 

\iPF = -, Pl = y^SN < f^SN , (11) 

^99 9 

when 

XpF = YpF = ^{l--fi) . (12) 
The pitchfork bifurcation sets in before the saddle-node bifurcation as can be readily seen from 

The upper branch of the homogeneous solution X^ given by ^ at the pitchfork bifurcation 
point fipF coincides with (jl2p . Hence as the Hopf bifurcation, the pitchfork bifurcation branches 
off the upper branch of the homogeneous solution. The pitchfork bifurcation is always subcritical 
because there are no solutions Xa possible for fj, > fipp- No bifurcation theory is needed to 
determine the subcritical character of the pitchfork bifurcation. 

The stability of the homogeneous solution X = Y = X^ = is determined by linearization. 
We study perturbations X = X^ + xexp At and Y = X^ + yexpXt, and obtain as a condition 
for nontrivial solutions x and y 

(A + 2gXh + /?7i) = ±/3e-^" . (13) 

The upper sign denotes an antisymmetric mode x = —y whereas the lower sign denotes a sym- 
metric mode X = y. Stationary bifurcations are characterized by A = 0, and in this case the 
symmetric mode coincides with the saddle-node bifurcation and the antisymmetric mode 
terminates at the pitchfork bifurcation (jl2p . 



As for the case of a single pulse in a ring, non-stationary Hopf bifurcations are possible li X = iu 
for wave trains. We obtain from ()13p 

uj = ^(3smiOT (14) 

and 

P 

Xh = —{± cos UJT - ji) . (15) 
25 

We consider only the symmetric case (the lower signs) which reproduces our results ([7]) and 
([8|) for the symmetry preserving Hopf bifurcation. The antisymmetric case does not allow for 
a single-valued positive uj. For ujt — > the onset of the Hopf bifurcation moves towards the 
saddle-node ([5]) and coalesces with it at /3r = 1 in a Bogdanov-Takens point as described in 
Section 12.11 For lot — > vr the Hopf bifurcation moves towards the pitchfork bifurcation with a 
limiting value of X^ = /3(1 — ji)/{2g) = Xpp in a codimension-2 bifurcation. At the point of 
coalescence the Hopf bifurcation has a period T = 2r which corresponds exactly to the inhomo- 
geneous pitchfork bifurcation with p = it whereby every second pulse dies. For values ujt € [0, vr) 
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Figure 3: Sketch of the bifurcation diagram for a wave train in a ring showing a stationary 
saddle- node bifurcation (SN) and a subcritical pitchfork bifurcation (PF). In between these two 
bifurcations is also a Hopf bifurcation (HH). 



the Hopf bifurcation always sets in after the pitchfork bifurcation. Hence for wave trains one 
can see only a Hopf bifurcation of the steady solution at the point where the Hopf bifurcation 
collides with the pitchfork bifurcation. 

This allows us to sketch the full bifurcation scenario for a wave train in a periodic ring as 
depicted in Figure [3l In the subsequent sections we study the bifurcations of the steady-state 
solution dl]). 

3 Bifurcation analysis of the Hopf bifurcation for a single pulse 

In this Section we study the direction and the stability of the Hopf bifurcation. In excitable 
media the Hopf bifurcation is a result of stronger and stronger coupling of the activator with 
its own inhibitor. Upon reducing the length of a ring for fixed excitability, or reducing the 
excitability for a fixed ring length, a single pulse will feel the tail of its own inhibitor created 
during its previous passage through the ring. In a wave train each pulse will feel the tail of the 
inhibitor of its neighbour in front. In the context of our normal form ([1]) the increasing coupling 
translates into an increase of /?r. Upon increasing /3r from /3t = 1 up to a critical value of 
/9r « 7.789 we have only one solution lo = loq oi the characteristic equation ([7]); see Fig. [6l 
We will study now the dynamics for this case. The case of arbitrary many solutions when one 
encounters a pseudo-continuum of frequencies will be discussed further down in Section HI 

The theory of bifurcation analysis for delay-differential equations is well developed [351 EH EH 
I38j . For example, in |38J a formula for the coefficients of the normal form for a Hopf bifurcation 
is given explicitly. However, we found that the theory of delay differential equations is not as 
well known amongst scientists as its age may suggest. We find it therefore instructive to perform 
the calculations explicitly and lead the reader through the calculations. 

To study the direction of the Hopf bifurcation we have to determine the sign of da/dfi at the 
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bifurcation point ^h- From ([6]) we infer 

d\ _ 2g dX 

dii~ 1 - Pre-^'^ djj, 



Using ([9]) we find from (jj]) 



and subsequently 



dX_ _ 1 

dl^ -coswr) ' 



d\ 2g 1 



,3 111- >°' 



which imphes that the stationary solution looses stability with increasing values of the bifurca- 
tion parameter ^. 

We now study the character of the Hopf bifurcation and derive the normal form for a Hopf 
bifurcation from ([1]). In order to do that we first transform the normal form ([1]) into standard 
form by subtracting the stationary solution ([4]) according to X = X + x where X = Xi. We 
obtain 

dtx = -2gXx - (3{x{t - r) + 7ix(t)) - gx^ , (17) 

with the stationary solution being now x{t) = 0. We will employ a center manifold reduction 
for this equation to describe the dynamics close to criticality. Center manifold theory is well- 
established for maps, ordinary differential equations and partial differential equations. However, 
although known for some time [35^ \3E\ l37] , it is not well known how to formulate an essentially 
infinite dimensional delay differential equation such as ([T|) into a form such that center manifold 
reduction can be applied. For ordinary differential equations, for example, the application of 
center manifold theory is a straight-forward expansion of the state vectors in critical eigenmodes. 
The problem for delay differential equations is their inherent infinite dimensional character. An 
initial condition x{9) = X(){9) for — r < 9 < is mapped onto a finite dimensional space; in 
the case of (jl7p onto a 1-dimensional space. Lacking uniqueness of solutions is one obstacle 
which prohibits a straightforward application of center manifold reduction. The trick out of 
this dilemma is to reformulate the problem as a mapping from an infinite-dimensional space of 
differentiable functions defined on the interval [— r, 0], which we denote as C = C[[— r, 0];M] (i.e. 
xo(^) € C), onto itself. This allows us to employ the well established and understood center 
manifold reduction for mappings. These ideas go back to Hale [MIES]. We found well written 
examples of center manifold reductions to be examined in [391 001 HI]. In essence, the history 
of a state vector x{t) G M is folded to a single element of an extended state space xt{9) £ C. In 
order to achieve this we define xt{9) G C as 

xt{9)=x{t + 9) for -T<9<0. 

The time-evolution for x{t) £ M (fT7|) needs to be reexpressed in terms of propagators and 
operators acting on elements of the extended state space xt{9) G C which can be done by writing 
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For — T < < we used the invariance condition dx{t + 6)/dt = dx{t + 6)/d0. For = we can 
split the operator T into a hnear part C and a nonhnear part A/" and reformulate the right-hand 
side of (fTTl) as 



where 



= £N+A^[xt] , (19) 



r-o 

C[xt]= I dewi{e)xt{e) with = -(25X + /37i)(^(^) -/3(5(0 + r) (20) 



and 

J\f[xt]= l\9^de2W2{di,92)xt{9i)xt{e2) with u;2(^i,^2) = -5<5(^i)5(^2), (21) 

where 5(0) denotes the Dirac 5-function. Once xt{0) is computed via solving (llSp . one may 
convert back to x{t) by 

.0 

x(t) = / de5{9)xt{0) . 
3.1 Linear eigenvalue problem 

We now linearize (jlSp around the stationary solution xt{9) = to obtain 



The linear eigenvalue problem ()22p can be solved using the ansatz 

^^(0) = e^*$(0) for -r<0<O. 
On the interval — r < < we obtain 

xm = , 

which is solved by 

^0) = e^^$(0) . (23) 



Plugging the solution (|23|) into (j22|) we can now evaluate the 9 = 0-part of (p2|) to obtain again 
the characteristic equation ([6]). We recall the transcendental characteristic equation as 

X + 2gX + P^i + Pe-^^ = . (24) 

Since in general the linear operator Al is not selfadjoint, we need to consider the corresponding 
adjoint eigenvalue problem on the dual extended state space Ct = Ct [ [0, r] ; M] . The dual problem 
is given by backward evolution for t < 0, i.e. xj(s) = X-t{—s) for < s < r. The adjoint problem 
can be formally written as 

j/ti^) = -Aims) . 
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To provide an explicit form of the dual operator Al we need to define an inner product. It 
turns out that the normal scalar product used for ordinary differential equations is not capable 
of respecting the memory effects of delay-differential equations. The following inner product is 
used 

(^t^$) = ^t(0)$(o) - y de j ds-^\s - d)wi{d)^{s) . (25) 
The adjoint operator is then explicitly given as 



A{ms) = \ -t4}'^ (26) 



where 

£t[^t] = / dswi{-s)il{s) . (27) 





The adjoint eigenvalue problem (|26|) is now solved using the ansatz 

^t(s) = e"^*^'^(s) for < s < r . 
On the interval < s < r we obtain 

which is solved by 

^'t(5) = g-As^t^O) _ (28) 

Plugging the solution (j28p into (j27p we can now evaluate the s = 0-part of (|26p to obtain again 
the characteristic equation Note that the two solutions ([25]) and ([25]) for the eigenvalue 

problem and its dual can be transformed into each other by simple time-reversal 9 —s. 
Since the transcendental equation has two solutions with vanishing real part a = 0, i.e. A = ztitu 
with oj given by ([7]) , we have two solutions of the linear eigenvalue problem (j22p and its associated 
adjoint problem (j26p . We denote them as $1^2 and 2 respectively. The bilinear form (j25p was 
constructed in order to assure biorthogonality of the eigenfunctions. Defining the eigenfunctions 
as 

$i(0)=e^'^^ ^2{9) = e-''^^ (29) 
^t(s) = ue-''^' ^{s) = v^e'^' (30) 

with the normalization constant 

we have (^'J, ^j) = 5ij with i,j = 1, 2 and 6ij being the Kronecker symbol. 
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Figure 4: (a): Stability diagram for the equilibrium solution We also impose (3t > 0. 
The region within the bold lines is stable. Crossing the upper boundary corresponds to a Hopf 
bifurcation, (b): Hopf line obtained by numerical simulations of the modified Barkley model 
(f70]) with a = 0.22 and = 0.1 and D = 3, jSSJ [Mj- (Note that each point corresponds to 
different values of e and L). The numerical results of the partial differential equations could be 
fitted to the normal form ([T]) to obtain the parameters /3o and r, see Ref. [28]. The continuous 
line is the same Hopfline as in (a). 



3.2 Center-manifold theory 

For the nonlinear theory we need properties for the linear operator A developed in [36\ I37| . 
We summarize properties of the transcendental characteristic equation ([6]): (i) .4 has a pure 
point spectrum, (ii) the real part of the eigenvalues is bounded from above, and (iii) defining 
a = T{2gX + (3^i} and b = (3t all eigenvalues of A have negative real part if and only if (1.) 
a > —1, (2.) o + 6 > and (3.) b < Csin^ — acosC where C is the root of C = — atanC with 
< C < TT for a 7^ and ^ = 7r/2 if a = 0. These conditions can be translated for our particular 
case (pH) using a = — (3t cos{ujt) . Condition (1) translates into /3rcos(a;r) < 1; condition (2) 
into cos(a;T) < 1 and condition (3) defines a parameterized stability boundary (3t < (/ sin{() 
and cos{ujt) > cos(C) with < C < ^r. This last condition defines the line in Pt-cos{ujt) space 
where the Hopf bifurcation occurs. In particular we have (3t > 1 and a coalescence with the 
saddle node at cos(a;r) = 1 with /?r = 1 at the Bogdanov-Takens point. In Figure H] we show 
the stability region. We note that care has to be taken in interpreting the diagram in terms of 
the parameter r because P = /3(t) according to ([2]). In the limit r — > oo we have /3r 0. Note 
that stable solutions may exist for /3r > 1. 

In [36^ [37] it is shown that under these circumstances one may perform center manifold 
reduction. In particular we can decompose xt{9) into slow modes associated with the eigenvalues 
A = ±iuj and fast modes which correspond to modes with negative real part of the eigenvalues. 
Center manifold theory says that the fast modes are slaved to the slow modes and can be 
expressed in terms of the slow modes. We therefore write 

xtie) = zmiie) + z''m2{0) + h{z, z'') , (32) 

where z{t) e C and its complex conjugate z*{t) are the time-dependent amplitudes of the slow 
modes $i,2(^) and h{z, z*) is the remaining fast component written as a function of the slow am- 
plitudes. The function h{z, z*) is called the center manifold. The expansion (132p is resemblant of 
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center-manifold theory for partial differential equations where would be the spatial coordinate, 
and the expansion would be an expansion of critical spatial eigenmodes. The connection between 
delay differential equations and partial differential equations will be explored further in Section^ 



We require the fast modes, and hence the center manifold, to lie in the spectral complement of 
the centre space spanned by ^{9); we therefore have the constraint 

= j = l,2. 

This implies for the slow amplitudes 

z{t) = {^>[{e),xt{e)) and z\t) = {m\{e),xt{e)) . (33) 

We use a near-identity transformation for the center manifold h and express it as a power series 
in z and z* . The center manifold is tangential to the manifold spanned by the slow modes which 
implies the ansatz 

h{z, z*) = ^ [h2o{9)z^ + 2hu{e)zz* + ho2{e)z*^) + 0{\zf) . (34) 

Since Xt{9) is real we have /io2(^) = ^20 (^)- Since the normal form for a Hopf bifurcation only 
involves cubic terms, we only need to consider quadratic terms in the equation for h (j37p . The 
cubic terms will then be generated via Af[xt] {0 = 0) in the equations for z and z* (see below, 
([35]) and ([36|)). 



We will derive now ordinary differential equations for the slow amplitudes z and z* which 
describe the dynamics on the slow manifold. The theory of center manifolds tells us that the 
full dynamics of (jl7p is well approximated by the slow dynamics |42) . The derivative of (j33p 
with respect to time t is given by 

i = {^{,Xt) 

= iu:{^\,xt) + {^,M[xt])\^^^ 

= iuJZ + ^l{0)M[xt]ie = 0) 

= iujz + i^M[xt]ie = 0) (35) 

z* = -iiuz* + i^*Af[xt]{e = 0) (36) 

h = xt-i<^i{e)-z*<^>2{e) 



= ALXt+M[xt] - iuJZ^i{e) + iuJZ*^2{0) - M[xt]{e = Q){u^l{e) + U*^2{0)} 

= ALh + N[xt] - {v^i{e) + u'^2{e)}M[xt]{e = o) , 



(37) 



where the dot denotes a time derivative. Note that A/'[x(] 7^ only for = which can be 
written as N[xt\{9 = 0) = N[z<^i + z*<^2 + K^, z*)]{9 = 0). Using ([Ml) we have therefore 



M[xt]{e = Q) = -9x1(9 = 0) 

= -5(z$i + z*$2 + /i)^|e=o 
= -giz"^ + z*"^ + 2\z\'^ 

+ /i2o(0)^^ + (/i2o(0) + 2hnm)\z\^z + (/io2(0) + 2hum)\z\^z' + /io2(0)z*^) 
+0{z\z*''). (38) 
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Using the definition of Al and Af[xt] we can evaluate the evolution equation ([18]) by differenti- 
ating the center manifold (|34p with respect to time, equate with (|37p . and obtain 

h = iuh2,{e)z^ - icvho2{e)z*' = -27^K"^]A/■[xi](0 = o) + 1 n ^ ^ ° (39) 



if6' = 
where 

n[z,z''-h] = -{2gX + P^i)h{0) - phi-r) +N'[xt]ie = 0) . (40) 
Here 7^ denotes the real part. 

Comparison of powers of z and z* yields differential equations for hij for the part with —t<6<0 
with an associated boundary value problem coming from a comparison of powers of z and z* 
from the 9 = part. We summarize 

h'20 = 2iLoh2o - AgTZiue''^^] (41) 
/i^2 = -2iioho2 - ^^gTlWe'"^^] (42) 
h[, = -457^K^], (43) 

and the boundary conditions are given by 

{iio - Pe-'^^)h2o{0) + (3h2o{-r) = 2g{2n[u] - 1) (44) 
i-iuj - pe''^^)ho2{0) + Pho2{-r) = 2g{2n[u] - 1) (45) 
i-iiv - (ie-'^^)hn (0) + (ihn (-r) = 2g{2n[u\ - 1) . (46) 

Note that the nonlinearity enters the differential equation in form of an inhomogeneity. The 
ordinary differential equations ()4ip -(H3]l can be solved by variations of constants 



h2o{0) = H2oe''^' -2i^[ve'^" + -v*e-''^'J (47) 

ho2{e) = Ho2e-^''^' + 2i^(^'e-'^' + ^^6'^'^ (48) 

hi^[e) = Hn + 2i^(ue^'^^ -iy*e-''^^^ . (49) 

The constants of integrations H2Q = and Hn can be determined using the boundary condi- 
tions (jMD-dM]). We obtain 

Note that Hn = —2g/{P — /3cos(wt)) = Hl^. By means of transformations [43j equation (135]) 
can be transformed into the standard form for a normal form for Hopf bifurcations 

i = iujz + c\z\'^z . (52) 



The quadratic terms appearing in (|35p can be eliminated by the near-identity transformation 
z ^ z + r]2oz'^ + + '?02-2*^ using 7/20 = igi^/^, "Hu = — 2zgz//u; and r/02 = —igi'/Suj. Note 
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Figure 5: Plot of the real part of the cubic coefficient 7^[c] (|53p as a function of ujt. To produce 
the plot we set g = 1 and r = 1. Both parameters are just coefficients multiplying c, so they do 
not change the sign of c. 



that at the Bogdanov-Takens point where a; = such an elimination of quadratic terms is 
not possible anymore. However this transformation generates further cubic terms in (j35p . All 
of these cubic terms except those proportional to \z\'^z may be eliminated by means of another 
transformation z — > z + h3{z, z*) where h^^z, z*) is a cubic polynomial. After the transformation 
to eliminate the quadratic terms the coefficient in front of the jzpz-term is found to be 



The stability and character of the Hopf bifurcation is determined by the sign of the realpart 
of c. Because of (fT6]) the Hopf bifurcation is supercritical provided TZ[c] < and subcritical 
provided Tl[c\ > 0. These criteria can be easily deduced by writing z = re*"^. Note that TZ[c] 
can be written as a function of g,T and lo only since (3 = to/ sin(u;T) at the Hopf bifurcation. 
Using algebraic software packages such as Maple, we can show that TZ[c\ > for all values of 
g, T and p. In Figure [5] we show the real part of c as a function of lot. This confirms that 
the Hopf bifurcation is indeed as conjectured in [28] subcritical. We have checked our result 
against numerical simulations of the full normal-form ([T]) and also using the software package 
DDE-BIFTOOL [3^. Again, the degeneracy at lot ^ is reflected in TZ[c\ by a singularity at 

UJT = 0. 

We will discuss the implications of this result in Section [5l 



c 



-5i^(/i2o(0) + 2/in(0)) 

+ {'^'no2gi^* + 2r]iigv* - r/ii(iLjr/2o + S^z^) - 2r/2o(2iwryii - gv) - 2r]Q2gv) 



-gu{h2oiO) + 2/in - — ig'— + ^^<?' — 
-gv{H2o + 2Hii) + ^ig"^— . 

6 LO 



(53) 
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Figure 6: Illustration of the solutions and number of solutions of the implicit equation (j54p for 
oj. Green curve: /?r = 0.1, no Hopf bifurcation; light blue curve: (3t = 5, Hopf bifurcation with 
one marginal mode; dark blue curve: fir = 143, Hopf bifurcation with finitely many marginal 
modes; pink curve: (3t = oo, Hopf bifurcation with infinitely many marginal modes. 

4 The limit of large delay times: The Hopf bifurcation for wave 
trains 

In the previous Section we have described the Hopf bifurcation for small ujt when there is only 
one marginal mode. This describes the behaviour of a single pulse on a ring. In this Section 
we will pursue the case of large delay times when a pseudo-continuum of critical Hopf modes 
occurs. We will derive a Ginzburg-Landau equation as an amplitude equation describing near- 
threshold behaviour of such a pseudo-continuum. The connection between amplitude equations 
and delay-differential equations has long been known [iSl HH HTJ HHl HHl |50] . The cross-over 
from a finite-dimensional center-manifold to an infinite-dimensional amplitude equation can be 
best viewed when looking at the Hopf condition ([7]) 

to = (3 sinujT . (54) 

For /3r > 7.789 there are at least two solutions of (j54p for uj. This equation has arbitrary many 
solutions LOk for /5r ^ oo and we obtain a pseudo-continuum in an interval with lower closed 
boundary at lot and upper boundary lot = vr. At the singular limit ujt = vr there are countably 
infinitely many eigenvalues cu^r = kn. An illustration is given in Fig. [6l Note that the upper 
boundary cut = vr corresponds to the coalescence of the Hopf bifurcation with the pitchfork 
bifurcation in the case of several pulses on a ring (see Section ^ . 

All these solutions are marginal and would have to be included in the ansatz (|32p . Note that 
for excitable media where /? = /9oexp(— «:r) (see ([2])) this limit cannot be achieved by simply 
letting T oo. The function /3t has a maximum at r = 1/k. So in order to have /3r — > oo 
one can either have /3o — oo which seems unphysical or k — > with r — s- oo to keep kt finite. 
Hence the limit /3r — > oo applies to media with a very slowly decaying inhibitor in a very large 
domain. Large domain instabilities are known from certain excitable reaction diffusion systems 
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in the context of autocatalytic oxidation of CO to CO2 on platinum |51[ 1521 [53]. 



The case ojt ~ tt is important for single pulses in a ring and for wave trains. Firstly, it describes 
the case for a single pulse when a continuum of modes becomes unstable to a Hopf bifurcation. 
But more importantly it describes the case of a wave train with distinct members when the 
pitchfork bifurcation coalesces with the Hopf bifurcation. The point of coalescence is at ^pf 
given by (jlip and with amplitude given by (jl2p which we recall 

(3 , 

^PF = 7^ 1 - 71 • 

At this point of coalescence the Hopf frequency is in resonance with the spatial instability in 
which every second pulse dies. At ^pp the two equations ([3]) describing the alternating modes 
in a wave train collapse to the single equation for one pulse ([l]) (see also Fig. [3]). This Section 
will investigate whether the coalescence of the subcritical pitchfork bifurcation with the Hopf 
bifurcation may produce stable oscillations. 

We will perform a multiple scale analysis of the normal form ([T|) along the lines of [48 j . We will 
obtain at third order an evolution equation for the amplitude as a solvability condition which 
describes the dynamics close to the Hopf bifurcation. We consider the case of large delay times 
r and introduce a small parameter e = 1/r. To capture the dynamics close to the point of 
coalescence we introduce a slow time scale 

s = et , 

and rewrite the normal form ([T|) in terms of the slow variable as 

edsX = -^i- gX^ - (3{-f + X(s - 1) + 71X) . (55) 
We expand the scalar field X{s) as 

X = xpp + exi + e^X2 + e^xs H . 

Using the generic scaling the bifurcation parameter can be written as 

/i = ixpp + e^A/x H . 

A Taylor expansion of (I54p around lot = tt yields at first order ujt = (5t{'k — lot) which for large 
T (small e) we may write as 

This suggests a multiple time scaling 

ds = dso + edsi + e^ds^ H . 

Close to the bifurcation point critical slowing down occurs which allows us to expand the delay 
term for large delays as 

X{s-1) = e-^'X{s) (57) 



.2/1 



''S2 



e'^'oX{s) . (58) 
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1 + e'^'o 



At lowest order, 0{1), we obtain the equation determining xpp. At the next order we obtain 

Cxi = , (59) 

with the Unear operator 

c = p 

Equation ([59]) is solved by 

xiisQ, si, S2) = z(si, S2)e^"^« + z{si, S2)e-*"^o , (60) 

with complex amplitude z and its complex conjugate z. Note that on the fast time scale t we 
would have xi{t) = zexp{iujt) + c.c. with ujt = vr, which, of course, is the Hopf mode at onset. 

At the next order, ©(e^), we obtain 

Cx2 = -Afi - gxl - ds^xi + pds^e'^'oxi . (61) 

The right-hand side involves terms proportional to exp(ibi7rso), which are resonant with the 
homogeneous solution of Cx2 = 0. We therefore impose the solvability condition 

dsoXi- pds.e-^'oxi = , 

which using (i59]l reads as 

ds^xi + ^dsgXi = . (62) 
In terms of the complex amplitude z using ([601) this reads as 

ds-t^z + iivrz = . (63) 

P 

This amounts to the time scale iujt « rdt ~ ds^ + edg^ = in — ien/f] = in^l — l/(/3r)) which 
corresponds to our scaling ([56j) at first order. Provided ([63]) is satisfied we can readily solve ([6T]l 
by solving for each appearing harmonic, and find 

X2 = [A;U + 2g|z|2+5z2g2i7rso^^^2g-2i^soj 

= -^[Af, + gxl] , (64) 

where we used ([59]) . 

At the next order, O(e^), we obtain the desired evolution equation as a solvability condition. At 
O(e^) we obtain 

Cx3 = -dsoX2 + Pds-.e-'^'o X2 - ds^xi - 2gxiX2 - '^(3ds^s^e~^''o xi + [5ds^e~'^''o xi 

-a 1 

= -dsoX2 + f3ds^e "0x2 - ds^xi -2gxiX2 + -Pds^siXi - Pds^xi . (65) 
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Again resonant terms proportional to exp(ibi7rso) are eliminated by imposing a solvability con- 
dition which upon using the expressions for X2 yields the desired amplitude equation 



(66) 



This is the well-studied real Ginzburg-Landau equation [54J- The time- like variable is the slow 
time scale S2 and the space-like variable the faster time scale sq which is 0{t). As in the finite 
dimensional case studied in Section [3] the Hopf bifurcation is clearly subcritical since the real 
part of the coefficient in front of the cubic term in (I66p is positive for all parameter values. Hence 
the coalescence of the Hopf bifurcation and the pitchfork bifurcation cannot lead to stable os- 
cillations. We have shown that wave trains also undergo unstable oscillations in the framework 
of the normal form ([1]). 

The usefulness of the spatio-temporal view point for delay differential equations as expressed 
here in the Ginzburg-Landau equation (i66l) has been pointed out @5l HHJ HU [50] . However the 
Ginzburg-Landau equation (I66p may be cast into a finite dimensional system which emphasizes 
the underlying multiple scale analysis. We start by rewriting (j66p as an equation for the complex 
amplitude z. One can explicitly express so-derivatives and obtain the following finite dimensional 
system 



The time-scaling on the left hand-side is as expected from our initial linearization and expansion 
of the frequency (I56p . We have in total 



which corresponds to (j56p . This illustrates the multiple-scale character of our analysis where 
the nonlinear term may be interpreted as a frequency correction |55j . The correction term to 
the linear term on the right-hand side of (j67p shows that the onset is retarded on the very slow 
time scale S2- 

In [21] a real Ginzburg-Landau equation was derived for paced excitable media with an addi- 
tional integral term modeling the pacing. It would be interesting to see whether the therein 
derived amplitude equation can be derived in a multiple scale analysis along the lines of this 
multiple scale analysis. 



5 Summary and Discussion 

We have explored the Hopf bifurcations of a single pulse and of a wave train in a ring of excitable 
medium. We have found that for the phenomenological normal form ([T|) the Hopf bifurcation 
for a single pulse on a ring and for a wave train on a ring is always subcritical independent on 
the equation parameters. 




(67) 
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Hopf bifurcations in excitable media had been previously studied. Besides numerical investi- 
gations of the Barkley model [33], the modified Barkley model [M], the Beeler- Renter model 
[Ml EH EH Hg EQl EH], the Noble-model [591 [HI |29] and the Karma-model [IZJ, where a Hopf 
bifurcation has been reported, there have been many theoretical attempts to quantify this bi- 
furcation for a single-pulse on a ring. Interest has risen recently in the Hopf bifurcation in the 
context of cardiac dynamics because it is believed to be a precursor of propagation failure of 
pulses on a ring. The Hopf bifurcation has been related to a phenomenon in cardiac excitable 
media which goes under the name of alternans. Alternans describe the scenario whereby action 
potential durations are alternating periodically between short and long periods. The interest 
in alternans has risen as they are believed to trigger spiral wave breakup in cardiac tissue and 
ventricular fibrillation fig [HI [13 [HI [1^ . 

Our results may shed a new light on what may be called alternans. The occurrence of alter- 
nans in clinical situations is often followed by spiral wave breakup and ventricular fibrillation 
|15 tllll[T6| ll7 t [T8j. The subcritical character of the Hopf bifurcation gives a simple and straight- 
forward explanation for this phenomenon. Moreover, if the system length L is slowly varied, 
long transients may be observed of apparently stable oscillations (see Figure [7| and Figure [8]). 
Depending on whether the system length is below or above the critical length Lh the oscillations 
will relax towards the homogeneous state or the instability will lead to wave breakup. However, 
even for the case of relaxation towards the stable homogeneous solution, these oscillations may 
lead to wave breakup upon further reduction of the system length, because of the subcritical 
character of the Hopf bifurcations. This illustrates the diagnostic importance of cardiac alter- 
nans. 



5.1 Limitations and range of validity of our results 

Strictly speaking, our result that the Hopf bifurcation is subcritical for the normal form ([T]) 
cannot be taken as a prove that alternans are unstable for all excitable media. The normal form 
([T|) is only valid for a certain class of excitable media. In particular it describes the situation in 
which an activator weakly interacts with the inhibitor of the preceding exponentially decaying 
inhibitor. Moreover, the normal form has only been phenomenologically derived in [28]. Of 
course, unless a rigorous derivation of the normal form ([1]) has been provided the results pre- 
sented here may serve as nothing more than a guidance in interpreting alternans in real cardiac 
systems or more complex ionic models of excitable media, and may alert scientists to check 
results on stability of oscillations more carefully. 

Several simplifications have been made to obtain the normal form ([!]) in [28]. For example, the 
time delay t = L/cq is treated as constant. This is obviously not correct for Hopf bifurcations. 
However, the inclusion of 71 (which is essential in the quantitative description of the Hopf bi- 
furcation) allows for velocity dependent effects. Guided by the success of the normal form to 
quantitatively describe a certain class of excitable media and by numerical experiments we are 
hopeful that our result may help interpreting experiments and numerical simulations. 

In Section [5.31 we will discuss a particular model for cardiac dynamics in which for certain pa- 
rameter values the assumptions for the derivation of our normal form are violated. For these 
parameter values stable oscillations may occur. However even for systems which are described 
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by the normal form ([T]) a word of caution is appropriate. If the oscihatory solutions bifurcating 
from the stationary solution are unstable as we have proven here, the unstable Hopf branch could 
in principle fold back and restabilize. Our analysis does not include such secondary bifurcations. 
Another scenario which we cannot exclude based on our analysis is that the unstable branch 
may be a basin of attraction for a stable oscillatory solution far away from the homogeneous 
solution. However, our numerical simulations do not hint towards such scenarios. 
From an observational perspective the relevance of the subcritical instability for spiral wave 
breakup is a matter of the time scale of the instability. The time scale associated with the 
subcritical Hopf bifurcation may be very long as seen in Fig. [9l This time scale becomes shorter 
the further the perturbation in the bifurcation parameter is from its value at the corresponding 
stable stationary pulse solution. In any case, if the parameter is kept fixed above the critical 
value, the instability will eventually develop unless the life time of a reentrant spiral is less than 
the time scale of the instability. For clinical applications one would need to estimate the time 
scale of a reentrant spiral and compare it with the time scale of the instability. Such estimates 
however are not meaningful for simple models such as the Barkley model. 

Our definition of alternans is restricted to non-paced pulses on a ring. If the excitable media is 
paced, the subcritical character of the Hopf bifurcation is not guaranteed anymore, and there 
is no a priori reason why stable alternans cannot occur. Indeed, in periodically stimulated 
excitable media stable alternans have been reported [60l EH [191 [SHI EH [23 [21] • A non-paced 
single pulse on a ring is a simple model for a reentrant spiral moving around an anatomical 
obstacle or around a region of partially or totally inexcitable tissue. As such it ignores the 
dynamics of the spiral away from the obstacle. An extension would be to look at a transversal 
one-dimensional slice through a spiral and consider wave trains and instabilities of such wave 
trains. 

5.2 Relation to the restitution condition 

Since the pioneering work [T2| alternans have been related to a period-doubling bifurcation. 
This work has rediscovered the results by [15], which had hardly been noticed by the scientific 
community until then. In there it was proposed that the bifurcation can be described by a 
one-dimensional return map relating the action potential duration {APD) to the previous re- 
covery time, or diastolic interval {DI), which is the time between the end of a pulse to the 
next excitation. A period-doubling bifurcation was found if the slope of the so called restitu- 
tion curve which relates the APD to the DI, exceeds one. A critical account on the predictive 
nature of the restitution curve for period-doubling bifurcations is given in [621 123j. In [29] the 
instability was analyzed by reducing the partial differential equation describing the excitable 
media to a discrete map via a reduction to a free-boundary problem. In |34j the Hopf bifur- 
cation could be described by means of a reduced set of ordinary-differential equations using a 
collective coordinate approach. In [Til EOl [SHI [26] the bifurcation was linked to an instability 
of a single integro-delay equation. The condition for instability given by this approach states - 
as in some previous studies involving one-dimensional return maps - that the slope of the resti- 
tution curve needs to be greater than one. However, as evidenced in experiments [631 El] and 
in theoretical studies [62l [231 ES [Ml 157] alternans do not necessarily occur when the slope of 
the restitution curve is greater than one. In our work we have a different criterion for alternans 
(which we interpret now as unstable periodic oscillations). Our condition for the occurrence of 
alternans, /3r > 1, does not involve the restitution curve but involves the coupling strength and 
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the wave length. Moreover, in Fig. Sl^b) we can see that for our normal form pulses can be stable 
for values of /3r ^ 1 in accordance with the above mentioned experiments and numerical studies. 

In the following we will show how our necessary condition for the onset of instability (3t > 1 
can be related to the restitution condition, that the onset of instability is given when the slope 
of the restitution curve exceeds 1. 



Close to the saddle node the Hopf frequency is ujt k. 0. We introduce a small parameter J <C 1 
and write close at the saddle node 

X = XsN + Sx , 

where Xsn is given by ([5]). The generic scaling close to the saddle node implies that we may 
write jj, = nsN + S'^^fJ-- Using the critical slowing down at the saddle node and the fact that 
UJT ~ we may approximate the normal form ([T|) to describe the temporal change of X at some 
time t and at some later time t + t. 

— — = -fJ-SN - <5^A^ - g{XsN + 6xnf - /3(7 + Xsn + IiXsn + Sxn-i + ii^Xn) ■ 

T 

Here x„ = and Xn+i = x{tn + t). Neglecting terms of and using the definition of 

the saddle node ([5]) we end up with 

Xn+l - (1 + /3t)x„ + l3TXn-l = . 

This equation has either the solution = 1 which corresponds to the stable steady solution 
described by Xi of ([1]), or 

Xn = (/3'r)"xo , 

which implies 

Xn = PrXn-l ■ (68) 

Close to the saddle node the amplitude of the activator correlates well with the APD, and we 
find that /3r > 1 is exactly the restitution condition whereby the slope of the restitution curve 
has to be larger than one. 



Our model contains the restitution condition as a limiting case when the Hopf bifurcation occurs 
close to the saddle node. However, as seen in Fig. H] /3r may be larger than one but still the 
system supports stable pulses. These corrections to the restitution conditions are captured by 
our model. Moreover, the normal form is able to determine the frequency at onset. 

We note that the parameter 71 does not enter the restitution condition; it is not needed for the 
existence of a Hopf bifurcation (cf. ([7|) and However, as pointed out in [28] quantitative 

agreement with numerical simulations is only given if 71 is included. In [S^ the inclusion of the 
7i-term takes into account the velocity dependent modifications of the bifurcation behaviour: 
large-amplitude pulses have a higher velocity than low- amplitude ones. A larger pulse will there- 
fore run further into the inhibitor generated by its predecessor. Velocity restitution curves have 
been studied in [65j to allow for a modification of the restitution condition derived in [30] for a 
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single pulse in a ring. The normal form incorporates naturally these velocity dependent terms. 

For a recent numerical study on the validity of the restitution condition the reader is referred to 
|67j . In this work the stability of certain excitable media is investigated by means of numerical 
continuation methods which allows a precise identification of the onset of oscillations. At the 
onset of alternans the restitution curve was determined. It was found that the restitution con- 
dition failed for three out of four cases for pulses in a one-dimensional ring. Our result suggests 
that the restitution condition may be a good indicator for the onset of alternans close to the 
saddle node. 



5.3 Numerical simulations 

In the context of alternans the Hopf bifurcation had been described as a supercritical bifurcation 
[T71 [29l [30j and not as we have found here as a subcritical bifurcation (although at the same 
time their occurrence had been related to wave breakup [29j). We therefore revisit some of the 
previous numerical studies. In [T7] the following two-variable model was proposed 



M' 

I n. \ 

A 



I n 



edtE = e''da:xE -E + 
dtn = e{E-l)-n, (69) 



Tp2 

(1 - tanh(£; - 3)) — 



as a model for action potential propagation in cardiac tissue. Here 9{x) is the Heaviside step 
function. This model incorporates essential features of electrophysiological cardiac models. For 
the parameters A = 1.5415, e = 0.009, M = 30 and ub = 0.525 a supercritical Hopf bifurcation 
was reported upon diminishing the system length L. We integrate this model using a pseu- 
dospectral Crank-Nicolson method where the nonlinearity is treated with an Adams-Bashforth 
scheme. We use a timestep of dt = 0.00001 and 4096 spatial grid points. A Hopf bifurcation 
occurs around L = 0.215. To approach the Hopf bifurcation we created a stable pulse for some 
large system length, and subsequently diminished the system length L. In Figure [7] we show that 
for these parameters the bifurcation is actually subcritical. The subcritical character has not 
been recognized before - probably because of insufficiently short integration times. For system 
length L just above the critical length the oscillations can appear stable for a very long time 
(see Figure [8]) before they settle down to the homogeneous solution. 

Indeed, as already stated in our paper [28], the number of oscillations may be rather large 
when the instability is weak. In Figure [9] we show such a case for the maximal amplitude of the 
activator u for the modified Barkley model 

dtu = DdxxU + u{l — u){u — Us — v) 

dtv = e{u-av) , (70) 

which is a reparameterized version of a model introduced by Barkley [32]. It is clearly seen that 
the oscillations can appear stable for a very long time and many oscillations (in this case more 
than 500 oscillations) which has lead scientists to the wrong conclusion that the Hopf bifurcation 
is supercritical. 

The normal forms ([T]) or ([3]) were derived for situations in which the activator weakly interacts 
with the tail of the preceding inhibitor which exponentially decays towards the homogeneous rest 
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Figure 7: Temporal behaviour of the maximal amplitude Emax of the activator E for model (169p 
just above the subcritical Hopf bifurcation. The parameters are A = 1.5415, e = 0.009, M = 30 
and ub = 0.525 and L = 0.215. The inlet shows the behaviour at L = 0.210. 




10 50 150 200 

Figure 8: Temporal behaviour of the maximal amplitude Emax of the activator E for model (j69p . 
The system length is just below the Hopf bifurcation with L = 0.22; the other parameters are as 
in Figure [71 (a): The oscillations appear to be stable over some time, (b): Same parameters as in 
(a) but longer integration time. The apparent stability has to be accounted for by insufficiently 
long integration times. The solution adjusts to the homogeneous solution. Note the long time 
scales which contain hundreds of oscillations. 
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Figure 9: Temporal behaviour of the maximal amplitude Umax of the activator u for model (j70p 
just above the subcritical Hopf bifurcation. The parameters are a = 0.22, Ug = 0.1, e = 0.03755 
and L = 246. The oscillations appear stable for a very long time but will eventually either damp 
out and attain a constant non-zero value in the case, when L is larger than the critical Lh at 
which the Hopf bifurcation occurs, or in the case L < Lh the pulse will collapse as depicted in 
Figure [7] confirming the subcritical character of the Hopf bifurcation. 



state. Then one can describe the influence of the tail of the preceding inhibitor as a perturbation 
to the generic saddle node of the isolated pulse. The models discussed so far all fall into this 
category. A different model was introduced by Echebarria and Karma in [21] which as we will 
see below for certain parameter regions does not fall into this class of model but supports stable 
oscillations. Originally the model was studied for a paced strand but recently has also been 
studied in a ring geometry [68j. It has been argued in ^68j that the stability of the spatially 
extended pulse is determined by the stability of a paced single cell. In the following we study 
numerically the Hopf bifurcation for this model in a ring geometry. This will illustrate the range 
of validity for our normal form and the conclusions which may be drawn with respect to the 
stability of cardiac alternans. The model consists of the standard cable equation 

dtV = Dd,,V-^^, (71) 



where lion models the membrane current and Cm is the capacity of the membrane. In [21] the 
following form for the membrane current was proposed 

^^ = -(s + {l-S)^)-^hS, (72) 

Cm To V ^cj Ta 



with a switch function 



2 V e 
The gate variable h evolves according to 

dh 1-S-h 



S=-(l + tanh(^^^— 5^) ) . (73) 



(74) 



dt Tm{l - S) + TpS 



The stable homogeneous rest state is at ^ = and h = 1; however for small a second sta- 
ble focus may arise. For details on the physiological interpretations of the model the reader is 
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Figure 10: Temporal behaviour of the maximal amplitude Vmax of the activator V for model 
dZU-dll]). The parameters are tq = 150, Ta = 26, = 60, Tp = 12, 14 = 0.1, D = 0.00025 
and e = 0.005. The main figure is obtained for L = 1.11 which is slightly above the subcritical 
Hopf bifurcation confirming the subcritical character of the Hopf bifurcation. The inset is for 
L = 1.1175 which is slightly below the bifurcation point. 



Figure 11: Space-time plot of stable oscillations occurring at Tq = 6 with L = 4.8. The other 
parameters are as in Fig. [TUJ Stable oscillations are found for a range of ring lengths L. 



referred to [21', 168] . For the numerical integration we use again a semi-implicit pseudospectral 
Crank-Nicolson method where the nonlinearity is treated with an Adams-Bashforth scheme. We 
use a timestep of dt = 0.01 and 1024 spatial grid points. In Fig. [10] we show an example for a 
subcritical Hopf bifurcation in this model consistent with our theory. However, for sufficiently 
small Ta a supercritical Hopf bifurcation arises upon decreasing the ring length L. In Fig. [TT] 
we present a space-time plot for such a situation of stable oscillations. Whereas the subcritical 
case is consistent with our theory we now have to understand why for small stable oscillations 
occur. In order to do so it is helpful to look at the spatial profiles of the activator and the in- 
hibitor close to the Hopf bifurcation which are presented in Fig. [T2l In the left figure we see the 
activator V and the inhibitor 1 — h for the case of a subcritical Hopf bifurcation as seen in Fig. 1101 
The figure is similar to Fig. [1] for the modified Barkley model. The activator weakly interacts 
with the exponentially decaying tail of the inhibitor it created during its previous revolution. In 
this parameter region our normal form is valid and correctly predicts a subcritical bifurcation. 
In the right figure of Fig. [12] the situation is depicted for the supercritical case seen in Fig. [TT] 
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Figure 12: Plot of the activator V (continuous line) and the inhibitor h (dashed line) for the 
system ()7ip -(|74 p . We plot here 1 — h rather than h to have the homogeneous rest state at 
u = and 1 — h = 0. Parameters are tq = 150, = 60, Tp = 12, Vc = 0.1, D = 0.00025 and 
e = 0.005. Left: The activator runs into the exponentially decaying tail of the inhibitor which 
decays towards the rest state 1 — h = 0. This is similar to the behaviour in Fig. [TJ Parameters 
are Ta = 26 with L = 1.11. This scenario is well described by the normal form. Right: The 
activator does not interact with the exponentially decaying tail corresponding to the rest state 
but rather with the metastable state defined by /i = 0. Parameters are = 6 with L = 4.8. 
This case cannot be captured by the normal form. 



Here the situation is very different. The inhibitor does not approach the homogeneous rest state 
l — h = but rather develops a metastable l — h = l plateau. This has two consequences; firstly, 
the solution is driven away from the homoclinic pulse solution around which the normal form is 
built, and secondly the interaction is not weak anymore. The reason for this different behaviour 
can be understood by looking at the nullclines of the homogeneous problem of (j7T|) - ([7H) . i.e. 
setting dx = 0. In Fig. [13] we show the nullclines for the two cases = 6 (supercritical) and 
Ta = 26 (subcritical) . Note that for Tq = 6 the only stable fix point is at ^ = and h = 1. The 
difference is that in the supercritical case the h = nullcline and the V = nullcline are very 
close to each other. This forces the trajectory to spend a long time on the h = 0-nullcline near 
h = (as seen in the plateau part of the spatial profile of 1 — /i = 1 in Fig. [TT]) . We call this 
state a metastable state. For decreasing ring length it dominates the profile of the inhibitor and 
does not allow the inhibitor to come close to the rest state h = 1. Therefore our normal form, 
which is formulated around the saddle-node of the pulse, breaks down. The solution is not close 
to the travelling pulse in phase space anymore and our local analysis around the saddle-node of 
the travelling wave cannot work anymore. However, we note that the system ()7ip -()74p is rather 
unusual with the two nullclines being parallel to each other with the possibility of a metastable 
state, resulting in rather particular dynamical behaviour. 

We may modify the model (j7ip - (|74p to break the degenerate situation in which the two 
nullclines of the inhibitor and activator run parallel to each other and subsequently may get too 
close to each other for certain parameters. We can destroy the existence of the metastable state 
for finite ring length by allowing the nullcline of the activator to bend away from the nullcline 
of the inhibitor if we, for example, consider the following modification of the membrane current 



with some sufficiently small ti . Then we are again in the situation where the rest state V = and 




(75) 
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Figure 13: Nullclines for the system (j7ip - (j74p . The continuous Unes denote the V = nullchnes 
and the dashed hues the h = nullchnes. Parameters are tq = 150, = 60, Tp = 12, Vc = 0.1 
and e = 0.005, and for all cases only one stable fix point exists at y = and h = 1. Left: The 
subcritical case with Tq = 26. Middle: The supercritical case with = 6. Note the closeness 
of the nullclines for large V. Right: Nullclines for the modification (j75p which breaks the near 
degeneracy of the nullclines observed in the middle figure. Here Tq = 6 and r; = 3. 



h = 1 dominates the dynamics upon decreasing the ring length L. The nullclines are shown in 
Fig. [T3j We confirmed that for r; = 3 the Hopf bifurcation is indeed subcritical, consistent with 
our theoretical result. We note that the actual value of r/ is not important for the existence of 
subcritical bifurcation but rather that a sufficiently small r; breaks the geometric structure of the 
degenerate nullclines and allows the activator nullcline to bend away from the inhibitor nullcline. 
The model (|7ip -(j74 p illustrates for which class of excitable media our normal form is applicable 
and for which systems we may draw conclusions on the stability of dynamical alternans in a ring. 
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